Intervaly spolehlivosti ve frekventistické interpretaci

frekventistické tvrzení má smysl, jen pokud lze experiment/událost opakovat!

95% horní mez $X_u$ NEZNAMENÁ, že hledaná (skutečná) hodnota bude s 5% pravděpodobností ležet nad ní
skutečná hodnota nějaké veličiny je dána (a neopakuje se)
tvrzení o horní mezi založené na měření (testovací statistice) $X$ znamená, že pokud je skutečná hodnota $X_u$, tak pravděpodobnost, že naměříme hodnotu $X$ nebo menší, je 5%

literatura: Barlow, sec. 7.2

Příklad

Z 20 pozorování nám 4 vyšla pozitivně. Jaký je 95% interval pro pravděpodobnost pozitivního výsledku?

Pravděpodobnost se řídí binomickým rozdělením. Příklad pro pravděpodobnost úspěchu p=0.3 dává 13% pro pravděpod. daného výsledku, necelých 11%, že vyjde méně než 4.



In [1]:

    
%pylab inline
from scipy import stats
vals=stats.binom(20,0.3).pmf(r_[0:15])
bar(r_[0:15]-0.4,vals)
title("binom. rozdeleni p=0.3")
vals[4],vals[:4].sum()
#stats.binom?









    



Populating the interactive namespace from numpy and matplotlib






    Out[1]:





(0.13042097437387012, 0.1070868045037307)



In [2]:

    
stats.norm.cdf(2)









    Out[2]:





0.97724986805182079

přehledné znázornění trojdim. grafu závislosti pravděpodobnosti na hodnotě $p$ a pozorované hodnotě $i$



In [24]:

    
x=r_[0.01:0.5:0.01]
for i in range(0,10):
    axhline(.1*i,color='gray',alpha=0.2)#,xmax=x[-1])
    plot(x,i*0.1+stats.binom(20,x).pmf(i))
axvline(.3)
xlim(x[[0,-1]])
xlabel("p")









    Out[24]:





<matplotlib.text.Text at 0x5309210>

pro (předpoklad) různé hodnoty $p$ lze spočítat pravděpodobnosti, že vyjde méně než daný počet (botlim) nebo více než daný počet (uplim) pozitivních výsledků...



In [11]:

    
botlim=lambda p,n:stats.binom(20,p).cdf(n-1)
uplim=lambda p,n:1-stats.binom(20,p).cdf(n)
x=r_[0.01:0.5:0.01]
plot(x,[uplim(a,4) for a in x],x,[botlim(a,4) for a in x])
legend(["uplim","botlim"])
xlabel("prst. uspechu")
ylabel("pravdep. obsah")









    Out[11]:





<matplotlib.text.Text at 0x4752f90>

Nyní zbývá najít řešení této tzv. Clopper-Pearson limity $$\sum_{r=m+1}^n P(r;p_+,n) = 0.975$$ $$\sum_{r=0}^{m-1} P(r;p_-,n) = 0.975$$ (tedy hodnoty $p_+, p_-$ jako horní a dolní mez konfidenčního intervalu)



In [12]:

    
from scipy import optimize
xbt=optimize.fsolve(lambda p:botlim(p,4)-0.975,0.1)[0]
xup=optimize.fsolve(lambda p:uplim(p,4)-0.975,0.4)[0]
xbt,xup









    Out[12]:





(0.057333997050032719, 0.43661400299666825)



In [14]:

    
#pozorovano 6 z 20
xbt6=optimize.fsolve(lambda p:botlim(p,6)-0.975,0.1)[0]
xup6=optimize.fsolve(lambda p:uplim(p,6)-0.975,0.4)[0]
xbt6,xup6









    Out[14]:





(0.11893159040572772, 0.54278918227628914)

Tahle interpretace je jediná smysluplná, pokud narazíme na fyzikální meze skutečných hodnot parametru (často jen kladné hodnoty, někdy omezené shora).



In [15]:

    
#pozorovano 6 z 20
xbt18=optimize.fsolve(lambda p:botlim(p,18)-0.975,0.6)[0]
xup18=optimize.fsolve(lambda p:uplim(p,18)-0.975,0.95)[0]
xbt18,xup18









    Out[15]:





(0.68301728598091771, 0.98765147282970489)

Cvičení

spočítat derivaci uplim(p) pro poissonovo rozdělení

Konstrukce pásu



In [26]:

    
n=100
cont=0.05
p=0.1
i=30
ptrue=r_[.05:1:0.05]
xlo=[stats.binom(n,p).ppf(cont) for p in ptrue]
xhi=[stats.binom(n,p).ppf(1-cont) for p in ptrue]
plot(xlo,ptrue)
plot(xhi,ptrue)
axvline(20)









    Out[26]:





<matplotlib.lines.Line2D at 0x9eff4ac>

Spojitá rozdělení

Pro normální rozdělení je situace jednodušší díky symetrii integrálu (graf z obr. 2 dává ve směru obou os v řezu gaussovky).

Podmínka $$\int_x^{\infty}{\frac{1}{\sigma\sqrt{2\pi}} e^{(x'-p_-)^2/2\sigma^2} dx'} =0.025$$ je ekvivalentní přímému výpočtu $$\int_{-\infty}^{p_-}{\frac{1}{\sigma\sqrt{2\pi}} e^{(x'-x)^2/2\sigma^2} dx'} = \int_{-\infty}^{p_-+x}{\frac{1}{\sigma\sqrt{2\pi}} e^{x''^2/2\sigma^2} dx''} = 0.025$$



In [32]:

    
xf=r_[:16:0.2]
qfun=stats.norm(20*0.3,sqrt(0.3*(1-0.3)*20))
bar(r_[0:15]-0.4,vals)
plot(xf,qfun.pdf(xf),'r')









    Out[32]:





[<matplotlib.lines.Line2D at 0x51e7ad0>]

aproximace úvodního problému rozdělením $N(\mu=n p,\sigma^2=n p(1-p))$



In [43]:

    
#qfun2=stats.norm(8*0.5,sqrt(0.8*(1-0.5)*8))
qfun2=stats.norm(4,sqrt(4*(1-4./20)))
qfun2.ppf(0.05)/20,qfun2.isf(0.05)/20









    Out[43]:





(0.052879819083977077, 0.34712018091602292)



In [40]:

    
# varianta prikladu 90% konfid. interval pri 4 uspesnych pokusech z 8
botlimq=lambda p,n,nn:stats.binom(nn,p).cdf(n-1)
uplimq=lambda p,n,nn:1-stats.binom(nn,p).cdf(n)
xbtq=optimize.fsolve(lambda p:botlimq(p,4,8)-0.95,0.1)[0]
xupq=optimize.fsolve(lambda p:uplimq(p,4,8)-0.95,0.4)[0]
xbtq,xupq









    Out[40]:





(0.1929029499941311, 0.80709705000589083)



In [42]:

    
qfun3=stats.norm(8*0.5,sqrt(0.8*(1-0.5)*8))
qfun3.ppf(0.1)/8,qfun3.isf(0.1)/8









    Out[42]:





(0.21343635827709961, 0.78656364172290039)